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Abstract 

Considered here are Boussinesq systems of equations of surface water wave theory over a variable bottom. A 
simplified such Boussinesq system is derived and solved numerically by the standard Galerkin-finite element 
method. We study by numerical means the generation of tsunami waves due to bottom deformation and we 
compare the results with analytical solutions of the linearized Euler equations. Moreover, we study tsunami 
wave propagation in the case of the Java 2006 event, comparing the results of the Boussinesq model with 
those produced by the finite difference code MOST, that solves the shallow water wave equations. 
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1. Introduction 

In recent years, there have been many theoretical and computational advances in the study 
of the full water-wave problem. However, there is still need for accurate, simpler mathematical 
models. Boussinesq systems are systems of partial differential equations that approximate the three- 
dimensional Euler equations that describe the irrotational, free surface flow of an incompressible, 
inviscid fluid. Because of their simplicity, Boussinesq systems have been used in the study of a 
variety of water wave phenomena in ports, channels, coastal areas, and in the open sea. They have 
also been used in studies of tsunami wave generation and propagation, cf., e.g., [TH], [15], [33] ■ 

In the case of horizontal bottom, these systems, derived in their general form in [3J, [5], may be 
written in nondimensional, unsealed variables as 

T) t + V • u + V • 77U + aAV • u - bArit = 0, ( , ( 

u t + V?7 + ±V|u| 2 + cAV?7 - dAu t = 0. 
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In these equations, the independent variable x = (x,y) represents the position, t is proportional to 
elapsed time, r\ = rj(x,t) is proportional to the deviation of the free surface from its rest position, 
while u = u(x, t) = (u±(x, t), 7/2 (x, t)) T is proportional to the horizontal velocity of the fluid at some 
height. The coefficients a, b, c, d are given by the formulas 

a = \ {d2 \ )v > 6 = \ {d2 - - v ^ c = V i - d= V i - ^ A* w 

where is, /u, are real constants and < 9 < I. (If we denote the nondimensional depth variable 
by z (with positive direction upwards), then the bottom of the channel lies at z = —1, while the 
horizontal velocity u is evaluated at the nondimensional height z = — 1 + 8(1 + rj(x, t)).) 

As it is explained in detail in [3J, the Boussinesq approximation on which (fT]) is based is valid 
when e := A/h -C 1, u := X/h ^> 1, and the Stokes number S = AX 2 //ig is of order 1; here A 
is the maximum free elevation above the undisturbed level of the fluid of depth ho and A a typical 
wavelength. Letting 5 = 1 and working in scaled, nondimensional variables, one may derive from 
the Euler equations, by appropriate expansion in powers of e, the scaled version of in the form 

r) t + V • u + e [V • ?yu + aAV • u - bArj t } = 0(e 2 ), ^ 
u t +Vr] + e [5VI11I 2 + cAV?/ - dAu ( ] = (9(e 2 ), 

from which ([1]) follows by rescaling and replacing the right-hand side by zero. 

It is worthwhile to mention that in [5] a new family of Boussinesq systems was derived. These are 
the fully symmetric Boussinesq systems of general form 



a,(uiit 2 )\ (4) 

vy ' \ + cAV? ? - dAu t = 0, 

where a, b, c, d are given by ([2]). These systems, like the usual Boussinesq systems ([lj, are formally 
0(e 2 ) approximations of the Euler equations when written in nondimensional, scaled variables in a 
form similar to that of ([3J), wherein their nonlinear and dispersive terms are multiplied by s. 

In addition, several types of Boussinesq systems with variable bottom have been derived, cf. e.g., 
[S], 0, S3], H3, EB, [12], US], [27]. The study of Boussinesq systems with variable bottom was 
initiated by Peregrine, [25] , who derived the system 

% + V-[(/i + £7 ? )u]=0 ! 

2 ( 5 ) 

u f + V?? + e(u • V)u - ct 2 |V(V • (hu t )) + ct 2 -^V(V • u t ) = 0{ea 2 , ct 4 ), 
where u denotes the depth-averaged velocity, i.e. 

u = / udz. 

h + er] J_ h 

Using ([5]) one may derive other Boussinesq systems. We mention, for example, Nwogu's system, [22], 

th+V- ((ev + h)u) + cr 2 V • (6 - I) /i 2 V(V • (hu)) + f^- - 6 + i) /i 3 V(V • u) = 0(ea 2 , ct 4 ), 
u t + V?? + e(u ■ V)u + a 2 \(9 - 1)W(V • (hu t )) + {9- 1) 2 ^V(V • u t )l = 0(e<r 2 , a 4 ), 

(6) 

where u denotes the horizontal velocity as in the case of the systems ([T]) . 
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In the paper at hand, we extend the Boussinesq models (0 and © in the case of a bottom that 
depends on x, y and t to classes of systems analogous to ([I) and ((4]). We also derive a simplified 
version of a BBM-type system appropriate for solving numerically realistic wave propagation prob- 
lems by the finite element method. Theoretical and numerical aspects of these systems in the case 
of a horizontal bottom, i.e., for the systems |T]) and (j4|), were studied recently in [5], [8], [10], pTj . 

We then apply the standard Galerkin-finite element method to the simplified Boussinesq model 
with homogeneous Dirichlet boundary conditions, and study the generation and propagation of 
tsunami waves. We compare the tsunami waves generated by this Boussinesq model and by the 
linearized Eulcr equations. Finally we study the propagation of a tsunami wave for a real tsunami 
event (that affected the island of Java in July 17, 2006), comparing the Boussinesq and the MOST 
models, [30]. The MOST model is an efficient numerical code solving the shallow water wave equa- 
tions combining a finite difference scheme and a splitting technique, cf. [30j . We do not study here 
the inundation caused by the tsunami, as the Boussinesq code is not yet equipped with a runup 
algorithm. 



2. Boussinesq systems over variable bottom 

2.1. Generalization of (Qp 

We denote by (i, y, z) a Cartesian coordinate system, with z measured upwards from the still 
water level. Consider a three-dimensional wave field with water-surface deviation propagating from 
its rest position, fj(x, y, t), at time t, over a variable bottom given by h{x, y, t) — D(x, y) + y, t). 
(We will assume that the variation of the time-dependent part of the bottom £ is of the same order 
of magnitude as the surface elevation fj.) The fluid velocity is denoted by u = (u, v, w) T . The Eulcr 
equations, which describe three dimensional wave propagation on the free surface, [32], are written 
in the form 

ur+(u- V)u+ivP= -gk, (7) 
P 

V-u = 0, for - h < z < fj(x,y,t), (8) 
V x u = 0, (9) 

where P is the pressure field, p is the density, g is the acceleration due to gravity, k = (0, 0, 1) T and 
V = (dx,dy,ds) T . The first equation expresses the conservation of momentum, whereas the other 
two equations express the conservation of mass and the irrotationality of the flow, respectively. The 
kinematic boundary conditions at the free surface and bottom can be expressed as 

f\i + u • V?7 = 0, for z = fj(x, y, t), (10) 

and 

h i + u-'Vh = Q, for z = -h(x,y,t), (11) 

respectively. The fluid is assumed to satisfy the dynamic boundary condition P = at the free 
surface z = fj(x,y,t). 

Consider a characteristic water depth ho, a typical wavelength Ao and a typical wave height ao, 
and the following nondimensionalization of the independent and dependent variables, cf. [25], [3], 
0, El], 

_ x_ _ y_ _ z_ _ co~ 

Ao Aq tlQ Ao 
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and 



-u, v 



h Q _ 

ti, w 



w, 77 = — , h=—, D = —, (=— , 

a c a /i ft a 



a c a c 

where Co = \f~gho~. Then, the governing equations for the fluid motion in nondimensional and scaled 
variables take the following form: 

1 



eu ( + e ((u • V)u + wu z 



p4 



VP = 0, 



1 



(12) 



ea 2 w t + e 2 a 2 (u ■ Vw) + ww z ) H ^P z = -1, (13) 

for —h < z < er\. Here u = (it, w) T and V = (d x , d y ) T , and the parameters e = cto/ho and cr = /iq/Aq 



are assumed to be small. The conservation of mass is formulated as 

V • u + w z = for — h < z < erj, 
while the irrotationality condition is expressed by the equations 

u y - v x = 0, 
u z - er 2 Vu? = 0, 

for —h < z < erj; the boundary conditions take the form 
rjt + e(u • V77) — w = on z = erj, 



(14) 

(15) 
(16) 

(17) 



Ct + u ■ V/i - 



on z 



(18) 



We note that h = D + e( and thus ht = 0(e). 

Now, following [25] (see also [13]), one may derive the generalization of Peregrine's equations with 
time-dependent variable bottom 



r\ t + V • [{h + £?y)u] + C* = 0, 

u ( + Vr; + e(u ■ V)u - <r 2 |V(V • (hu t )) + a 2 ^V(V • u t ) - <7 2 f VC ti 



0(£ct 2 ,c7 4 ). 



(19) 



Following the methodology of [52] (see also [12] and [3j), consider the horizontal velocity of the 
fluid u e at some height z = —h + 6(erj + h), with < 6 < 1. Then, one may derive from (IT2|) ~ ([T8]) . 
using appropriate expansions, the generalization of ^ given by 



r)t + V • (hu e ) + eV • (rju e ) + cr 2 V • a/i 2 V(V • (hu e )) + bh 3 V(V • u< 

+a 2 aV ■ {h 2 V( t ) + Ct = 0(ea 2 , a 4 ), 
u| + V?? + e(u e • V)u e + a 2 [cW(V • (/m?)) + Jfr 2 V(V • uf)l + a 2 chVC, tt = 0(ea 2 ,a 4 ), 



(20) 



where a = 6< - \, b = \ [(9 - l) 2 - |] , 5 = - 1 and J = ±(0 - l) 2 . (We note that the system (gDJ) 
is Nwogu's system with z Q = (0 — l)h in the notation of [22] and u = u e + 0(er 2 ).) 

We observe that due to the irrotationality condition (fT5|) . (u 9 • V)u e = iV|u e | 2 + 0(a 2 ), and 
V(V ■ u e ) = Au 9 + 0(cr 2 ). Moreover, there holds that hAu 6 = V(V • {hu e )) - V(V/i • u e ) - V/iV • 
u e + 0(a 2 ), and V(V ■ (/iuf)) = V(Vh ■ u 9 t ) + VhV ■ u 9 t + hAuf + 0(a 2 ). So, the system after 
dropping the superscript 6, may be written in the form 
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r\ t + V • (hu) + eV • (rju) + cj 2 V • |(a + b)h 2 W(W ■ (hu)) - Wi 2 [V(V/i ■ u) + V/iV • u]| 

+a 2 5V-(/ l 2 VC t ) + C t = 0(ea 2 , ( T 4 ), 
u t + V77 + e±V|u| 2 + (j 2 jc7i[V(V/i • u t ) + V/iV • u t ] + (c + d)/i 2 Au t } 
+a 2 chVC,tt = 0(e<r 2 ,a 4 ). 
From these equations we observe that 

^ = -V-(/ui)-Ct + 0(£,<7 2 ), and u t = -V V + 0{e, a 2 ). (22) 
Let fti/Gl, Splitting the dispersive terms using ([22]) (in the spirit of 3]) as follows 

V(V • (/m)) - /xV(V • {hu)) - (1 - M)V(r? t + 0) + 0(e, a 2 ), 
Au t = vAu t - (1 - v)AVri + 0(e, a 2 ), 
and ignoring the high order terms, we may write the system (|21[) in the form 

r) t + V • ((ft + er?)u) + cr 2 V • {Ah 2 [V(Vh ■ u) + V/iV • u] + aft 2 V(V • (hu)) - bh 2 Vn t } 

+a 2 AV ■ (h 2 VC t ) + (t = 0, 
u t + V?/ + e±V|u| 2 + a 2 {Bh[V{\7h ■ Vrf) + VftAr?] + ch 2 V{Ar]) - dh 2 Au t } - o- 2 BhV( tt = 0, 

(24) 

and in dimensional form 

r] t + V • ((h + r))u) + V ■ {Ah 2 [V{Vh • u) + VftV • u] + aft 2 V(V ■ (hu)) - Wi 2 Vr7 t } 

+iV-(/i 2 VC t ) + Ct=0, (25) 

u t + gVr? + ±V|u| 2 + {Bgh[V(Vh • V??) + V/zAr?] + cg/i 2 V(Ari) - d/i 2 Au t } - BftVCt* - 0, 
where 

A=^-(9-l)% 5 = 1-0, A = //a-(l- M )6, 

c = ±(i-0>, d=i(i-tf»)(i-„). 

We note that the parameters a, b, c, d are those of the class of Boussinesq systems derived in [3]; 
if we take the depth h as constant, the systems (|24|) reduce to the analogous systems of [3]. 

2.2. Generalization of 

Following the same procedure as in [5] one may generalize (j4|) to a class of Boussinesq systems 
similar to (|25|) . Specifically, considering the shallow water wave equations 

rit + V • ((h + eri)u) + ( t = 0, 

( 26 ) 

ut . + .9V7/+ §V|u| 2 = 0. 

and using the nonlinear change of variables hu — u(h + £77), the fact that u = u + O(e), and that 
/it = e£tj one may derive in dimensional variables the system 



T) t + V • (hu) + ±V • (ur?) + V • {A/i 2 [V(V/i • u) + VW • u] + ah 2 V(V ■ (hu)) - bh 2 Vrj t } 

+iv • (^ 2 vc t ) + Ct = o, 

u * + + JhdWv + sV|u| 2 + 27TUV • (hu) + ±u( t 

+ {Bg[V(Vh ■ Vt?) + VhArj] + cgKV(Ari) - dhAu t } - cr 2 BS7( tt = 0, 

where A, A, B, a, b, c, d as before. If we neglect the dispersive terms of (|27|) . then the resulting system 
conserves the function E(t) = L 2 gij 2 + h\u\ 2 , in the sense that E(t) = £7(0). Moreover, if we choose 
constant depth h, we recover the fully symmetric Boussinesq systems derived in [5]. As Peregrine 
pointed out in [35], h and its derivatives must be of O(l) to ensure the validity of all the above 
models. For other symmetric Boussinesq systems over variable bottom we refer to [B]. 

Remark: For specific examples of systems produced by specific choices of the parameters fi, 
is, 8 we refer to [3J, 0], [2], [5]- We mention that other triplets of parameters may be chosen so 
that the dispersion relation of the Boussinesq system approximates better the dispersion relation of 
the full water wave problem. Such a choice is for example, // = 0, v = (25 — Vl605)/49 and 9 2 — 
(80-\/l605)/105, i.e. a = 0, b ^ 0.0235121, c S -0.0952381, d = 0.4050592, which gives a (2,4)-Pade 
approximant of the full water wave problem. Another choice is /1 = —0.3672365, v = —0.3301459 
and 6 2 = 0.3906251 (i.e., a = -0.0105198, b = 0.0391657, c = -0.1005913, d = 0.4052788), which 
gives a (4,4)-Pade approximant of the full water wave problem. 



3. The simplified system and the numerical method 



In computations we have used a simplified Boussinesq system of BBM type. Consider the system 
(g5D with 9 2 = 2/3 and fi = v = 0, i.e., A = i= y / |-|, B = 1- J\, a = c = Q,b = d= 1/6. 
Then assuming that the bottom h is such that the derivatives of order greater one are of 0(e) 
and omitting terms of 0(e 2 ,£<J 2 ), we consider the following initial-boundary- value problem for the 
system (|25p for 1 in a plane bounded domain and t > 0: 

■q t + V ■((£> + C + »?)u) + 2ADVD ■ V(V • {Dm)) - 6V • {D 2 V m ) + AV ■ (D 2 V( t ) + Ct = 0, 

(28) 

u t + 5V77 + iV|u| 2 + BgD[(VD ■ V)Vr/) + VDA77] - dD 2 Au t - BDVC, tt = 0, 

n{ x iV^) = Vo(x,y), u(x,y,Q) =u (x,y), v(x, y, 0) = v (x, y), (x,y)eCl, 
with Dirichlet boundary conditions r) = u = v = 0ort dQ, for t > 0. 

We solved the above initial-boundary-value problem using the standard Galerkin-finite element 
method with continuous PI elements on a triangulation of fi. For the time stepping we used an 
explicit Runge-Kutta method of order 2 with a uniform timestep. This numerical scheme was an- 
alyzed and used in [TU] and [TT], in the case of the BBM-BBM Boussinesq system with constant 
depth. The main difference in the case of a general variable bottom is that the matrices that the 
semidiscretization of the u-equations yields are not symmetric, and thus for the numerical solution 
of the linear systems for u and v we use the Generalized Minimal Residual Method (GMRES) with 
appropriate ILUT preconditioner, cf. |26j . For the bottom h we use the interpolant in the finite 
element space Vh, while we approximate the initial data using an appropriate elliptic projection into 
Vh- For example, the initial condition for the u component of the solution is the function uo,h G Vh 
which satisfies 

-4("o,h,x) = -4(1*0, X), for all X ^ V h , 
where A is the bilinear form defined by the relation 
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A(u, v) = (u, v) + d(D 2 Vu, Vu) + d(X7D 2 Vu, v), for all u,v &Hq, 

where Hq is the usual Sobolev space Wq ,2 {VL). 

From the previous comments it is obvious that we will have to impose restrictions on the smooth- 
ness of the bottom topography to ensure the well-posedness of (|28| and the stability of the numerical 
method. For the triangulation of the domain f2 we usually use the "Triangle" software, [28] . 



4. Tsunami generation phase 



Many tsunami waves are caused by a bottom dislocation near a rupture due to an earthquake. In 
that case the bottom deformation may be approximated by Okada's formulas, [23], [23]. We briefly 
describe the vertical component of Okada's formulas in the case of the so called dip-slip dislocation, 
which is the most important component in the case of tsunami generation. 

In this case, consider a rectangular fault of width W and length L positioned near the z < 
axis and at a depth d below the free surface z = 0. The vector D represents the slip on the fault. 
The dip angle 5 and the angle 4> between the the fault plane and the slip vector, describe a general 
dislocation of the fault, where the vertical component of the displacement vector, (which we will 
denote here by 0{x,y)) 1 is given by the following formulas in Chinnery's notation, /(£, rj)\\ : = 
f(y + ^p)-f(y + i,p-W)-f(y-±,p) + f(y-±p-W),cf.mi,m- 



0{x,y) = ~ 

Z7T 



— ; + sin d arctan 1 sin o cos o 

R(R + qR 



(29) 



where U 
R 2 



|D|sin< 



e 
i = 



, p = y cos 5 + d sin 5, q = y sin 5 — d cos 5, y — i] cos S + q sin 5, d = rj sin S — q cos 5, 
= £ 2 + y 2 + d 2 , X 2 = £ 2 + q 2 . When cos£ ^ 0, / is given by the formula 



fi 



arctan 



r/(X + qcos6) + X(R + X)smd 
£{R + X)cos5 ' 



[i + A cos (5 
and, when cos 6 = 0, by 
j /i £ sin S 
M + Ai? + J' 

Here (i, A are the Lame constants given by the formulas /i = E/2(l + v), X — Ev/(1 + v){\ — 2v). 
The parameter E is the Young's modulus and v is the Poisson's ratio; both are considered to be 
known constants. In the sequel, it is supposed that Okada's vertical deformation component 0(x 1 y) 
has been translated to the appropriate location. 

In [14] and [17] , some mechanisms of the dynamics of tsunami generation are described. We review 
here two cases. A common practice in the literature is to use as initial condition for the free surface 
elevation Okada's solution, while the initial velocity profile is considered to be zero. This is referred 
to as passive generation of a tsunami and is described by the initial conditions 

r,(x,y,0) = 0{x,y), u(z,y,0) = 0. (30) 

Another way to model the generation of a tsunami is by active generation. In this case we consider 
zero initial conditions for both the surface elevation and velocity field, and assume that the bottom 
is changing in time. This case may be described by considering the bottom motion formula 

h(x,y,t) = D(x,y) + 0(x,y)F(t), (31) 

where T is a function of time t. In this paper we will consider three cases used in [T3], [T3J, [17] . 
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T t {t) :=H{t), 

T e (t) := 1 - exp- Kt , k > (32) 

T c {t) := H(t - t ) + |[1 - cos(7rt/t )]H(<o - t). 

In the numerical experiments that follow we use k — 2 and to — 2sec in the above. 

In each of the above cases the analytical solution of the linearized Euler equations has been 
computed in [T3] and is given by the formulas (see also [T7]): 

m{x,y,t) = ^J R 2 C( t'oih(mfe + r coButdkdt, 

Ve(x,y,t) - j^yz j R2 cosh(mh) ^ ^2+^2 ) ami, 

r) c (x,y,t) = (2^/ R2 c 'osh(mh) (sin ujt ~ cos jt + H(t - t )[cos uj(t - t ) + cos -/t]) dkdi, 

where 7 2 = j-, m = \Jk 2 + £ 2 , u> — \/ gm tanhm/i. The symbol C represents the Fourier transform 
of ((x,y,t), i.e. 

C(fc,^)= / ax,y,t)e-^ x +^dxdy. 
Jm 2 

In the numerical experiments we evaluated numerically the above formulas by approximating the 
usual Fourier transform by the discrete Fourier transform, using the FFT. 

In the sequel, we present a comparison between the numerical solution of the Boussinesq model 
(f2"5| and the analytical solution of the Euler equations (|33[) to study the generation of a tsunami. 
For the deformation of the bottom we used the formulas (j2"9j) - (f32)) in the square [—2,2] x [—2,2] in 
geographical coordinates (i.e. longitude-latitude coordinates in degrees). More precisely, we consid- 
ered the bottom motion given by h(x, y, t) = D(x, y) + 0(x, y)T(t), where T is one of the functions 
(f3"2")l . and D(x,y) = D , with D = 500, 1000 and 3000m. To compute the vertical displacement 
function 0(x,y) from (|2^)) we use the set of parameters shown in Table [TJ Some comparison results 
appear in figures [THSJ These figures show the surface elevation produced by the two models as a 
function of spatial variable x when y = in three cases. The horizontal scales in these figures are in 
geographical coordinates (degrees). The spherical shape of the earth was taken into account, even 
if it does not play significant role because of the small spatial scale of the experiments, (near the 
equator one degree is approximately equal to 111km.). 



Table 1 

Parameters used for the first experiment 



Parameter 


Value 


Parameter 


Value 


Dip angle <5(deg) 


7.0 


Fault length L (km) 


40 


Rake angle </>(deg) 


67 


Fault width W (km) 


20 


Strike angle (deg) 


90.0 


Fault depth d (km) 


10 


Slip amount, |_D|(m) 


2 


Young modulus E (GPa) 


9.5 


Longitude (deg) 





Poisson ration v 


0.27 


Latitude (deg) 





Acceleration of gravity g (m/sec 2 ) 


9.81 



For the finite element code we used 90970 elements, while for the computation of the FFT we 
used a rectangular grid of 57600 squares. (In these experiments we chose the parameters such that 
both models arc valid. Recall that the Boussinesq model is valid when the Stokes number S is 0(1). 
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t = 500sec t = lOOOsec 

Fig. 1. Do = 500m, case T e , k = 2. Boussinesq — , Euler T £ , Euler T c — ■ — , Euler Ti 





V. 







t = lOsec 




t = 200scc 



t = 400sec t = 600sec 

Fig. 2. Dq = 3000m, case T c , to = 2. Boussinesq — , Euler T e , Euler T c — ■ — , Euler JF^ 
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t = lOsec t = lOOscc 




t = 200sec t = 500sec 

Fig. 3. Dq = 3000m, case jFj. Boussinesq — , Euler T e , Euler T c — ■ — , Euler T% ■ ■ ■ 

In practice this means that 0.5 < S < 35, cf. pQ. For the choice of the time scales in formulas (|32l) . 

cf. HZ].) 

We observe that the solution of the Boussinesq model is close to the analytical solution of the 
linearized Euler equations (|33|) . The best agreement is apparently achieved in the cases where we 
used T c and Ti. For a more thorough study of tsunami generation we refer to [12] . 



5. A case of tsunami propagation: Comparison between MOST and the Boussinesq 
code 

In this section we present the results of a simulation of the propagation of the tsunami wave 
that affected the island of Java in July 17, 2006, using the Boussinesq model ([28} and the MOST 
code. MOST is an efficient numerical model solving the nonlinear shallow water wave equations in a 
characteristic form, using a splitting technique, coupled with an explicit second-order in space and 
first-order in time finite-difference scheme, [25], [5UJ, [3Tj . 

The seismological data that we used for the tsunami generation were provided by the National 
Oceanic and Atmospheric Administration - NOAA 1 . The fault's length and width were 100km and 
50km, respectively, the strike angle was taken equal to 199deg, while the epicenter was placed at 
(107.18,-9.56) in geographical coordinates. (Poisson's ratio was v — 0.23 in this case.) For the 
bathymetric data we used the GEBCO one-minute grid provided by the British Oceanographic 
Data Centre - BODC 2 . For the definition of the coastlines which represent part of the boundary 
of the domain f2, we used the Global Self-Consistent Hierarchical High-Resolution Shoreline data, 
GSHHS 3 . 



1 http://www.pmel.noaa.gov 

2 http:/ /www. bode. ac.uk/projects/gebco. html 

3 http:/ /www. ngdc.noaa.gov/mgg/shorelines/gshhs. html 
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In the case of the Boussinesq model, the bathymetric data were interpolated and smoothed appro- 
priately to ensure the stability of the numerical method and the validity of the Boussinesq system. 
However, because of the large variations of the bottom (see Figure [4]), shorter waves were gen- 
erated, especially around Christmas Island (southwest of Java) and around the undersea canyon 
near the earthquake's epicenter. In addition, due to some incompatibility between the shoreline and 
the bathymetric data we ignored riffs and islands that were not included in the shoreline data by 
changing the depth close to the shoreline to 100m. In the case of MOST we used the original bathy- 
metric data. (We note that MOST defines the coastlines by an internal procedure when a depth 
less than 10m is detected.) In both models we used the passive approach for the generation of the 
tsunami, i.e., initial conditions (l30|) . The time step for both models was taken equal to lsec. For the 
Boussinesq code we used 85863 elements while for MOST we used a uniform grid of width equal to 
O.Oldeg. (We ran the same experiment using the Boussinesq code with 80628 and 181480 elements 
with no significant differences.) 




103 104 105 106 107 108 109 110 



Fig. 4. Bathymetry of the sea around the island of Java used by the Boussinesq model, and position of the wave 
gauges and the epicenter of the earthquake. 0: wave gauge, □: the epicenter. (All distances in degrees.) 

In addition to surface elevation contour plots, we measured the variation of r/ as a function of 
time at eight wave 'gauges' placed at the positions represented by <0> in Figured! Specifically, gauges 
were placed at the points (i) (107.5,-8.5), (ii) (106.5,-8), (iii) (108,-8.5), (iv) (107.3,-8.8), (v) 
(106.9,-10.6), (vi) (107.7,-11), (vii) (105.9,-10.35), (viii) (104.9,-10.8). The results are shown 
in Figure [5J 

In Figure[5]we present the evolution of the initial data and the propagation of the ensuing tsunami 
in a series of surface elevation contour plots. We observe that the tsunami waves produced by the 
Boussinesq code and MOST model have similar shapes except in the dispersive tail (see also Figure 
[5]). The rightmost images of Figure [5] represent the maximum of r\ at (a;, y) up to t — 1750sec. (The 
color scale in the last two graphs of Figure[5]is between the values —0.3 and 0.5 but we checked that 
the solution, especially near the coastline, exceeds the value of 0.9m.) From these graphs one may 
observe that in the case of the Boussinesq model the front of the tsunami wave is more narrowly 
directed than the front produced by MOST. In Figure [6] we observe that the amplitudes of the 
solution of the Boussinesq and the MOST models agree quite well at all the gauges. We conclude 
that in this experiment the results of propagation of the tsunami wave produced by the MOST 
model and by the Boussinesq model are quite similar. For more information about the specific event 
we refer to [IB] . 
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Boussinesq model 
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MOST model 
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Fig. 5. T] at t = 0, 500, 1000, 1500sec, and maximum of r\ up to t = 1750sec. (All distances in degrees. Wave heights in meters). 



Gauge (i) 



Gauge (ii) 




Gauge (vii) Gauge (viii) 

Fig. 6. Tj-component of the solution at the eight wave gauges shown in Figure [4] as a function of t. Boussinesq 
MOST - -. 
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